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ABSTRACT 

A recent ALMA image revealed several concentric gaps in the protoplanetary disk surrounding the 
young star HL Tau. We consider the hypothesis that these gaps are carved by planets, and present a 
general framework for understanding the dynamical stability of such systems over typical disk lifetimes, 
providing estimates for the maximum planetary masses. We collect these easily evaluated constraints 
into a workflow that can help guide the design and interpretation of new observational campaigns 
and numerical simulations of gap opening in such systems. We argue that the locations of resonances 
should be significantly shifted in massive disks like HL Tau, and that theoretical uncertainties in the 
exact offset, together with observational errors, imply a large uncertainty in the dynamical state and 
stability in such disks. This presents an important barrier to using systems like HL Tau as a proxy 
for the initial conditions following planet formation. An important observational avenue to breaking 
this degeneracy is to search for eccentric gaps, which could implicate resonantly interacting planets. 
Unfortunately, massive disks like HL Tau should induce swift pericenter precession that would smear 
out any such eccentric features of planetary origin. This motivates pushing toward more typical, 
less massive disks. For a nominal non-resonant model of the HL Tau system with five planets, we 
find a maximum mass for the outer three bodies of approximately 2 Neptune masses. In a resonant 
configuration, these planets can reach at least the mass of Saturn. The inner two planets’ masses are 
unconstrained by dynamical stability arguments. 

Subject headings: celestial mechanics— planets and satellites: dynamical evolution and stability 


1. INTRODUCTION 

Given its wavelength coverage, sensitivity and resolu¬ 
tion, the Atacama Large Millimeter/Submillimeter Ar¬ 
ray (ALMA) represents a unique observatory for study¬ 
ing not only the extent and masses of protoplanetary 
disks, but also their structure. The latter provides the 
possibility of detecting forming planets as they gravi¬ 
tationally sculpt the surrounding disk. From October 
24-31 2014, as part of its science verification process, 
the ALMA team pointed 25-30 antennas at the young 
star HL Tau (J2000 04:31:38.45 +18:13:59.0), measuring 
the continuum emission in Band 6 (211-275 GHz), with 
baselines from 15 nr to 15 km 5 . As of the time of this 
writing, no scientific publication has appeared analysing 
these data; however, the ALMA project has publicly re¬ 
leased an image of the system (see Fig. 1) that reveals 
several concentric dust gaps, perhaps revealing forming 
planets massive enough to sculpt the nearby disk mate¬ 
rial. 

This striking discovery during the science verification 
phase suggests that ALMA may find similar features in 
several other young systems. In this paper, we there- 
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fore generally investigate the constraints one can extract 
on potential planetary masses from the requirement that 
observed systems be dynamically stable. We then apply 
these constraints to the HL Tau system using approxi¬ 
mate gap separations from the publicly released image. 
As argued in Appendix A, our analysis should be largely 
insensitive to an improved determination of the orbital 
radii. 

We begin in Sec. 2 with a list of assumptions that we 
make in our analysis. Section 3 then investigates stabil¬ 
ity in non-resonant systems and applies our thresholds 
to numerical simulations of the HL Tau system, finding 
that resonances may be important in shaping the sys¬ 
tem’s stability. Section 4 discusses the role of resonances, 
and provides simple criteria for evaluating whether res¬ 
onances must be considered for stability. Section 5 out¬ 
lines a procedure for numerically investigating resonant 
systems, applying it to HL Tau. Section 6.1 provides a 
summarized workflow for how one might analyze a po¬ 
tential planetary system embedded in a gas disk, making 
reference to more detailed discussion in the main text. 
Our conclusions for the HL Tau system can be found in 
Sec. 6.2. 


2. ASSUMPTIONS 

By specializing to protoplanetary disks observed by 
ALMA, we can restrict the relevant range of scales, sim¬ 
plifying our analysis. For any pair of putative planets 
(throughout the paper, subscripts 1 and 2 to the inner 
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and outer body, respectively) we assume 


Aa 1 

(1) 

— > 10 - 1 , 

ai 

p i ,2 > 10 yrs, 

(2) 

\< 10 - 2 , 

M* 

( 3 ) 

t < lMyr 

( 4 ) 


where a, P, and M represent the semimajor axis, orbital 
period, and mass, respectively; A a is the orbital sepa¬ 
ration, the * subscript refers to the central star, and t 
is the age of the system. Since the maximum resolution 
of ALMA is > 10 milliarcseconds, and the closest star 
forming regions are > 100 pc away, the smallest struc¬ 
tures ALMA can resolve are a few AU in size (« 5 AU 
in the HL Tau image). If we assume that the maximum 
orbital radius of planets is ~ 100 AU, Aa/ai > 10 -1 . In 
the case of HL Tau, the minimum separation between the 
putative planets is Aa/ai ss 0.2. If we assume most ob¬ 
served stars will be less massive than the Sun (M 0 ), then 
planets in resolvable gaps (beyond a few AU from the 
star) will have orbital periods > 10 yrs. Assumption (3) 
restricts us roughly to planetary masses. At about this 
threshold, corrections to the below results from terms 
that are higher order in the planet-star mass ratio be¬ 
come important. 

We additionally assume that the observed system’s 
age is at most a few Myr old, as this is the median 
timescale on which gas disks disperse (e.g. Haisch et al. 
2001; Cloutier et al. 2014). We note that all our max¬ 
imum masses are still valid over longer timescales, they 
just become less relevant. On Gyr timescales, a range 
of more subtle secular and chaotic effects can destabilize 
progressively lower-mass systems (see, e.g. Lithwick and 
Wu 2011; Wu and Lithwick 2011; Batygin et al. 2015). 
The short timescales in question (< 10 5 orbits) allow us 
to provide particularly simple estimates. 

Finally, by integrating the current system forward in 
time, we are implicitly assuming that the current state 
is representative of its past (or at least migrated into the 
current state over a timescale comparable to the system’s 
age). In particular, we assume planets did not recently 
jump to their current positions. This is reasonable in the 
HL Tau system where the orbital periods of the outer 
planets are ~ 10 3 yrs, since simulations of gap opening 
in the dust component of gaseous disks (Fouchet et al. 
2010) show that the timescale for the planet to clear a 
gap is 10 — 100 orbits. 

3. INSTABILITIES IN NON-RESONANT SYSTEMS 


3.1. Two planets 


Despite the three-body problem notoriously not pos¬ 
sessing enough constants of motion to be integrable, 
Marchal and Bozis (1982) and Milani and Nobili (1983) 
showed that, for a large enough combination of total en¬ 
ergy and angular momentum, the motion is constrained 
in such a way that close encounters can never occur (i.e., 
the system is Hill stable). For initially circular and copla- 
nar orbits, Gladman (1993) translated this stability cri¬ 
terion to a separation of A a > 3.461?#, where Rh is 
their mutual Hill radius 


Rh = «i 


\ 1/3 

M\ + M2 \ 

3M* J 


( 5 ) 


Planets separated by less than 3.5 Rh destabilize on 
the timescale of conjunctions t con j (where planets arrive 
at the same longitude and close encounters occur). For 
ALMA systems, this corresponds to t con j < 10 orbital 
periods of the innermost detected planet 6 —much shorter 
than the disk’s lifetime. 

In most two-planet cases, separation greater than 3.5 
Rh (and thus avoidance of close encounters), is a suf¬ 
ficient condition for stability (Barnes and Greenberg 
2006). It is nevertheless possible to more slowly diffuse 
through chaotic regions of phase space and lose planets 
(e.g., Veras et al. 2013); however, this is unlikely to play 
an important role for systems probed by ALMA. The 
chaotic region where this is possible corresponds closely 
to the Hill-unstable region, and only extends above 3.5 
Rh for planets with M < 10(M*/M 0 ) Earth masses 
(M 0 ) (Deck et al. 2013). These masses are smaller than 
those generally found capable of producing perturbations 
on the dust and gas distributions in protoplanetary disks 
that would be observable by ALMA, e.g., Dong et al. 
2014. Additionally, given our range of orbital periods 
(Eq. 2), Veras et al. (2013) find that such evolution is 
generally slower than the lifetimes of typical protoplane¬ 
tary disks (Haisch et al. 2001). 

Therefore, for (non-resonant) two-planet systems ob¬ 
served by ALMA, we conclude that Aa > 3.5 Rh fur¬ 
nishes a robust criterion for stability. For an observed 
separation, and assuming equal-mass planets, one can 
invert the equation to obtain the maximum mass for sta¬ 
bility M crit , 

m<m - = 1 ! (|)(tt) (6) 

For unequal mass planets, the condition is simply that 
the sum of the masses be less than twice M cr j t . 


In this section, we begin by considering the simpler 
case of a disk-free system, assuming that resonances be¬ 
tween planets are not modifying the dynamics. In sub¬ 
sequent sections, we then sequentially incorporate ad¬ 
ditional layers of complexity and investigate how they 
modify our simple estimates. 

Much work has been done on the minimum separation 
between planetary orbits that is required for stability, 
as a function of mass. We first discuss the case with 
two planets orbiting a central star (i.e., the three-body 
problem), which exhibits qualitatively different behavior 
from the many-body problem. 


3.2. Three or more planets 

The addition of a third planet expands the degrees of 
freedom available to the system, such that conservation 
of the total angular momentum and energy can no longer 
guarantee the avoidance of close encounters for initial 
separations beyond 3.5 Rh- A detailed theoretical un¬ 
derstanding of stability in systems with three or more 

6 tconj = 27r/An where An is the difference in the bodies’ mean 
motions n = 2tt/P. For An/m 1, Kepler’s 3rd law yields 
tconj ~ (Aa/ai) -1 orbital periods. For ALMA systems, Aa/ai > 
0 . 1 . 
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planets is currently lacking; however, numerical simula¬ 
tions furnish useful results. 

Chambers et al. (1996) found that systems of three 
terrestrial-mass planets, spaced by the same number of 
Hill radii from one another, all go unstable, at least for 
spacings Aa/Rn < 10. However, they found that the 
timescale to instability scales exponentially with A a/R h- 
Similar behavior was also observed for giant planets by 
Marzari and Weidenschilling (2002) and Chatterjee et al. 
(2008), though the scalings quantitatively differ. Nev¬ 
ertheless, the above studies agree that for separations 
< 4,Rh, instability sets in after less than ~ 10 3 orbits. 
This translates to the conservative stability threshold' 

M s = 7 8 (^)(^r) M ®- (7 » 

Despite this estimate not being analytic like Eq. 6, Eq. 7 
should represent a maximum mass in the sense that it 
applies to the idealized case of planets equally separated 
in A a/a. One would expect that as one separated the 
third planet from the other two, one would recover the 

2- planet condition (probably not formally, but certainly 
over the relevant timescales). Marzari (2014) recently nu¬ 
merically studied such equal-mass but unequally spaced 

3- planet systems. His results suggest that when the sep¬ 
aration (in A a/a) of the third planet is « 3 — 4 times 
larger than that between the tighter pair, one roughly 
recovers the two-planet result. Thus, in a general sys¬ 
tem, if one identifies the minimum A a/a, Eqs. 6 and 7 
bracket M CI in the limits of an isolated pair of nearby 
planets (a hierarchical system) and of a chain of equally 
spaced (in A a/a) planets, respectively. 

3.3. The effects of eccentricity damping 

Protoplanetary disks will additionally perturb plane¬ 
tary orbits (Papaloizou and Terquem 2006), and this will 
modify the results of the previous two sections. These in¬ 
teractions are complicated—they depend on the detailed 
balance between the gravitational effects of gas at res¬ 
onant locations with the planet, while the gas densities 
at these locations depend on the disk’s thermal struc¬ 
ture and are in turn modified by the resonant torques. 
As a result, the planetary mass required to open a gap 
in the disk (e.g. Bryden et al. 2000; Crida et al. 2006), 
and the associated timescales on which planetary or¬ 
bital elements evolve, are uncertain. Rather than ar¬ 
gue for specific scenarios or timescales, we choose simple 
parametrizations for how the orbital elements will evolve, 
and consider a wide range of timescales. One can then 
use our results to set dynamical constraints by selecting 
the appropriate timescale under a particular scenario, 
e.g., a planet that opens a deep gap, a planet that opens 
a gap in the dust, but not the gas, etc. 

We begin by considering the effects of eccentricity 
damping by the disk (Goldreich and Tremaine 1980; 
Artymowicz 1992, 1993). 8 Following Papaloizou and Lar¬ 
wood (2000) and Teyssandier and Terquem (2014), we 

7 Our definition of the Hill radius differs from that used in the 
studies above. We chose the definition given in Eq. 5 because it’s 
the one appropriate in the analytic limit of Eq. 6. In any case, the 
discrepancies are small, and always act in the direction of making 
our limit (Eq. 7) even more conservative. 

8 We do not consider the possibility of eccentricity excitation by 


implement eccentricity damping by including an addi¬ 
tional force in our N-body simulations 9 , 

F e = -—rf, (8) 

7~e 

where F e is the eccentricity-damping force vector acting 
on a planet, r is the radial component of the velocity, 
and r is the unit vector pointing outward in the radial 
direction. This yields an exponential damping of eccen¬ 
tricities on a timescale r e . This assumes that eccentricity 
is damped at constant angular momentum. 

We now qualitatively assess the effect of eccentricity 
damping on the stability of planetary systems. As argued 
above, if a pair of planets is separated by less than 3.5 Rh, 
the system will destabilize on a timescale < 10 orbital pe¬ 
riods. Therefore, unless the eccentricity damps on orbital 
timescales (and objects so tightly coupled to the gas seem 
unlikely to sculpt the disk into an observable feature), 
Eq. 6 furnishes a stringent upper limit on the planetary 
masses, independent of the disk’s properties and/or the 
depth of the gap the planet opens. For separations be¬ 
yond 3.5 Rh, the gas-free orbital evolution is more grad¬ 
ual; we would expect that if the eccentricity damping 
timescale is shorter than the instability timescale (which 
scales exponentially with A a/Rn, Sec. 3.2), the disk will 
suppress instability. 

Thus, we conclude that even in the case of three planets 
equally spaced in A a/a (Sec. 3.2), one should approach 
the two-planet result (Eq. 6) in the limit of strong eccen¬ 
tricity damping. We show this in our simulations of the 
HL Tau system in the next section. 

3.4. Application to HL Tau 

We now apply the above constraints to the HL Tau 
system and compare our estimates to numerical integra¬ 
tions. Not knowing the stretch in the publicly released 
image, we chose to focus on the most significant gaps, 
labeled in Fig. 1. We discuss in Sec. 3.5 how additional 
planets change our results. We describe how we deter¬ 
mined the orbital radii corresponding to each gap, and 
provide a table of values in Appendix A. Additionally, 
to limit the parameter space, we assume that putative 
planets all have equal mass. 

Our nominal radii of 13.6, 33.3, 65.1, 77.3 and 93.0 AU 
(table 1) correspond to relative spacings between adja¬ 
cent gaps (a,; + i —afj/ai of 1.4, 1.0, 0.19, and 0.20. Since 
the smallest separation is shared by two adjacent pairs, 
we are in the limit of Eq. 7. Additionally, as discussed 
above, since the inner two putative planets are > 5 times 
more separated than the outer three, the former are un¬ 
likely to affect the dynamics. Indeed, though the analy¬ 
ses described below include all five planets, at our grid 
resolution, we obtained indistinguishable results for the 
system’s stability with only the three outer planets. In all 
our integrations, we adopt M* = 0.55M 0 , estimated by 
Beckwith et al. (1990) using protostellar evolution tracks 
on luminosity-temperature plots. 

the disk (Goldreich and Sari 2003; Tsang et al. 2014), since (at least 
in closely packed systems) this would rapidly lead to instability. 

9 A library for parametrized migration, damping of eccentrici¬ 
ties and inclinations, and disk-induced precession can be found at 
https://github.com/dtainayo/rebound/tree/migration (see the 
readme on that page). It can be used to call various REBOUND 
integrators from within python. 
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Fig. 1. Publicly released continuum 233 GHz image of HL 
Tau taken during ALMA’s science verification process. The gaps 
we identified and used in our simulations (see Appendix A) are 
displayed. Blue lines aid centering the ellipses along the major and 
minor axes. The dashed ellipse shows our assumed outer limit of 
the disc, which we used to set our distance scale (Appendix A). 
Taken from http://www.eso.org/public/news/esol436/. 

We therefore estimate from Eq. 7 that the planets must 
have M < 30M®. Additionally, in the limit of strong ec¬ 
centricity damping (with a timescale somewhat longer 
than the conjunction time of ss (0.2) -1 orbits ~ a few 
thousand yrs), we expect that the slower 3-body insta¬ 
bility will be suppressed, and that we should recover the 
2-planet limit preventing close encounters (Eq. 6). This 
corresponds to M < 45M®. 

For our numerical simulations, we used the high-order 
integrator IAS15 (Rein and Spiegel 2015), which is part 
of the REBOUND package (Rein and Liu 2012). We 
define instability as when a planet’s radial excursions 
become greater than 5 AU, either through semimajor 
axis variations, or through radial variations (ae) induced 
by an excited orbital eccentricity. Like previous studies 
(e.g., Gladman 1993), we find our results are not sensitive 
to exactly how we define our criterion. 

In Fig. 2 we show the time to instability as the masses 
of the planets and the eccentricity damping timescale 
are varied. Each grid point corresponds to 24 numeri¬ 
cal integrations lasting 10 6 yrs (the upper-limit for the 
star’s age), where the particles were initialized at ran¬ 
dom longitudes on effectively circular and coplanar or¬ 
bits (e = 1CT 8 ,* = 10 -8 rad ~ 10 -6 deg, where i is the 
orbital inclination). The color in each grid point repre¬ 
sents the median of the time to instability (or 10 6 yrs for 
runs that remained stable) across the 24 simulations. We 
also tried runs with 96 simulations per grid point and ob¬ 
tained visually indistinguishable results. Along the top 
of Fig. 2, we use the mass to express the separation be¬ 
tween the closest planets in units of mutual Hill radii (see 
Sec. 3). 

We see that the system is somewhat more unstable 
than our estimates would suggest. While this is not 
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Fig. 2.— Time to instability for our nominal 5-planet case, 
as a function of planet mass, and eccentricity damping timescale. 
The top axis shows the corresponding separation between the clos¬ 
est planets in units of Hill radii. Each grid point is the median 
timescale to instability from 24 numerical integrations (see text). 
The vertical blue lines labeled 2 PL and 3+ PL correspond to 
the predictions expected in the strong eccentricity-damping and 
dissipation-free limits at 3.5 and 4 Rh, respectively (see Eqs. 6 
and 7). 

necessarily cause for alarm (Eq. 7 provides a maximum 
mass— lower masses may still be unstable), the main rea¬ 
son for the discrepancy is that the outer two planets are 
near the 4:3 resonance (their period ratios are ss 1.32). In 
the next section, we discuss how such resonances mod¬ 
ify the simple picture painted above. Nevertheless, we 
see the effect expected from the previous section that for 
decreasing r e , one approaches the 3.5 Rh threshold. 

Because of the sensitive dependence of stability on sep¬ 
aration, one simple possibility that would render the sys¬ 
tem substantially more stable is if not all gaps corre¬ 
spond to planets. For example, gaps 3 and 4 may not 
represent two separate planets, but rather a single one, 
with co-orbital material on horseshoe orbits in between. 
We therefore also consider a system with four planets at 
13.6, 33.3, 71.2 (the average of gaps 3 and 4) and 93.0 
AU (Fig. 3). This corresponds to A a/a = 1.4,1.1 and 
0.31. Given the single close pair of planets, and inner 
two planets at much wider separation, we would expect 
the stability limit to lie close to 195 M® (Eq. 6); how¬ 
ever, the two outer planets’ period ratio is now ss 1.49, 
near the 3:2 resonance, which pushes the stability limit 
to lower masses (because of our e « 0 initial conditions— 
see Sec. 4). 


3.5. Dependence on initial conditions 

Before investigating the role of resonances, we briefly 
justify our choices of initial conditions, and discuss how 
they affect our results. We first address our effectively co¬ 
planar geometry. One would only expect inclinations to 
matter on these timescales if they are capable of mitigat¬ 
ing close encounters. This should roughly occur when the 
maximum height reached on an inclined orbit a sin / be¬ 
comes comparable to the separation between planets Aa, 
since beyond this inclination, strong encounters would 
only occur when conjunctions happened to line up with 
one of the nodes where the two planets’ orbital planes 
cross one another. However, we argued (Eq. 1) that res¬ 
olution is likely to limit ALMA to finding separations 
A a/a > 0.1. This corresponds to i > 10 deg, which is 
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Fig. 3.— Time to instability for a 4-planet case where gaps 3 
and 4 are the result of a single planet with co-orbital material in 
between. Otherwise same as Fig. 2. 


high if the disk effectively damps inclinations (e.g., Pa- 
paloizou and Larwood 2000). In any case, our brief ex¬ 
perimentations with such large inclinations suggest that 
planets avoid ejections through the additional constraint 
that close encounters must occur near the nodes; how¬ 
ever, they tend to evolve onto crossing orbits that would 
not produce the observed concentric rings. 

Redoing the simulations in Fig. 2, but now drawing in¬ 
clinations from Rayleigh distributions with i = 10~ 4 and 
10 -2 rad (« 0.006 and 0.6 deg), both with and without 
inclination damping, we find visually indistinguishable 
results from Fig. 2. This motivated our choice for the 
nominal scenario presented in Sec. 3.4 with essentially 
coplanar initial conditions, ignoring inclination damp¬ 
ing. We briefly note that, while Barnes et al. (2015) 
found that planets in resonances with small inclinations 
could chaotically evolve onto extremely eccentric orbits 
that would certainly destabilize a tightly packed system, 
this occurs on much longer timescales than we are con¬ 
sidering. 

As for the eccentricities, one might expect that the ini¬ 
tial values are not significant as long as they are much 
smaller than the eccentricity required for crossing orbits 
Across ~ A a/a, which is approximately 0.2 in our HL 
Tau simulations. Indeed, if we redo the analysis in Fig. 2 
(which has initially circular orbits), drawing the various 
eccentricities from Rayleigh distributions 10 with scale pa¬ 
rameters of 10~ 4 and 10~ 2 (and the longitudes and peri- 
centers randomized), we obtain visually indistinguishable 
results from Fig. 2. However, as we discuss in the next 
section, this is likely not the appropriate way to choose 
initial conditions when planets are near resonance. 

We also point out that our assumed mass of 0.55M 0 
(Beckwith et al. 1990) is uncertain. Deviations from this 
value would affect our result; however, solving for a sep¬ 
aration of a certain number of Hill radii (Eqs. 6 and 7) 
yields the mass ratio between star and planet. Thus, 
our values (and errors) scale directly with the star mass. 
The secondary effect that the orbital periods will change 
is not as significant, given the exponential rise in insta¬ 
bility timescale with separation. 

10 with the Rayleigh probability density function / defined 
through f{e: cr) = xexp[—x 2 /(2 <t 2 )]/<t 2 , where a is the scale pa¬ 
rameter. 


Finally, we comment on the inclusion of additional 
planets to our nominal case. There may be two faint 
gaps between gaps 2 and 3 in our nominal model (Fig. 1), 
as well as a possible gap near the outer edge. If we sim¬ 
ulate as above a system of these 8 planets (with orbital 
radii of 43.4, 53.4 and 107.2), we obtain qualitatively 
similar results to Fig. 2, though shifted over one bin to 
lower masses at our level of resolution. The small discrep¬ 
ancy reflects the fact that the additional planets do not 
introduce much smaller A a/a values (though the out¬ 
ermost two planets have a somewhat smaller A a/a of 
0.15). However, one important difference is that, with 
the increased number of planets, it becomes harder to 
find orbital radii within the errors that are not close to 
some first-order mean motion resonance (MMR) (while 
in the next section we easily find such configurations for 
the nominal five-planet case). Resonances are therefore 
important in an eight-planet scenario (see next section). 
We note that this does not necessarily imply an eight- 
planet scenario requires lower masses—the reduction in 
stability is a product of our choice of initial conditions 
(Sec. 4). 

4. RESONANT EFFECTS 

Resonances can be perplexing in that they sometimes 
provide great stability (e.g., the 2:3 resonance between 
Pluto and Neptune that has prevented collisions over 
several Gyr despite their crossing orbits), while in oth¬ 
ers cases they are destabilizing (e.g., the Kirkwood gaps 
at resonant locations with Jupiter in the asteroid belt). 
We now briefly review the relevant dynamics in order to 
identify the relevant regimes and provide a lens through 
which to interpret numerical investigations. 

4.1. Background 

On the short timescales we are interested in, systems 
destabilize through close encounters at conjunctions. Be¬ 
cause resonances modify both the planets’ orbital eccen¬ 
tricities, as well as where the conjunctions occur, they 
can have a strong effect. As we discuss below, whether 
these changes help or hinder stability depends on ini¬ 
tial conditions, so these must be judiciously chosen. For 
studies of stability near the 3:1 and 2:1 resonances on 
longer timescales without a gas disk, see Marzari et al. 
(2005, 2006). 

It is well known that, near an isolated MMR between 
two planets on nearly circular and co-planar orbits, one 
can obtain an approximate one-degree-of-freedom Hamil¬ 
tonian that can be analytically integrated (e.g., Message 
1966; Peale 1986). Sacrificing precision for simple es¬ 
timates, we draw intuition from the case of a massless 
particle perturbed by a massive planet on an exterior, 
coplanar and circular orbit. 11 We additionally special¬ 
ize to the case of first-order mean motion resonances 
(with period ratios of the form j + 1 : j). which dom¬ 
inate when the eccentricities are small (see Murray and 
Dermott 1999, Chap. 8). 

11 For a particularly clear (but involved!) treatment of the more 
general case with two massive bodies on elliptical orbits, see Deck 
et al. (2013). They show that one obtains the same dynamics as 
in the test-particle case for first-order MMRs, when expressed in 
terms of a generalized eccentricity that is a weighted sum of the two 
planets’ orbital eccentricities. Our simple treatment thus provides 
rough but qualitatively accurate estimates. 
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In this case, one obtains a solution in terms of the ec¬ 
centricity of the test particle e± (the outer orbit is circular 
and fixed) and the resonant angle <f> = jX\ — (j+ 1)A2—roi, 
where A and w correspond to the particles’ mean lon¬ 
gitudes and their longitudes of pericenter, respectively. 
Physically, one can think of the resonant angle as trac¬ 
ing where conjunctions will occur, as when the planets 
line up (Ai = A 2 ),</> = A — i.e., <j> is the longitude 

of conjunction, measured from pericenter. This is an 
important angle since, for example, if conjunctions al¬ 
ways happen at pericenter, this maximizes the distance 
between particles; as a result, </> = 0 is the stable config¬ 
uration for first-order MMRs. Exact resonance occurs at 
semimajor axes that render <fi stationary, i.e., 


<j> = jn 1 - (j + 1 ) 77-2 - mi = 0, (9) 

where n represents the mean motion 2n/P. This can be 
trivially rearranged to give 


P 2 _ j + l _ 

Pi j jn 2 ' 


( 10 ) 


Generally, w\ will precess slowly compared to the plan¬ 
ets’ orbital rates, so exact resonance approximately cor¬ 
responds to the intuitive separation where P 2 /P 1 = 
(j + 1 )/j; however, we will see below that this is not 
always the case. 

A useful way to visualize the dynamics is as a polar 
plot of ei and </>; see Fig. 4. A point on this plot spec¬ 
ifies the state of the system at some time, with the ec¬ 
centricity given by the point’s radial distance from the 
origin (the dashed gray circles of constant eccentricity 
are labeled at e = 0.01,0.02,0.03), and (j) by the angle 
from the positive x axis. For a given planet separation 
that yields periods approximately satisfying Eq. 10, one 
finds that along <f> = 0, there is a stable equilibrium for 
a particular value of the eccentricity ej, referred to as 
the forced eccentricity. For initial conditions slightly dis¬ 
placed from the equilibrium point, the system traces out 
a closed contour around this fixed point, corresponding 
to oscillations in time for both the eccentricity and <fi. 

On short timescales, resonances can either help or hurt 
stability, depending on initial conditions. Superposed on 
the smooth evolution shown in Fig. 4 are discrete kicks 
at conjunctions. If the system starts sufficiently far from 
the equilibrium, the potentially large rises in eccentric¬ 
ity, together with the variations in conjunction longitude, 
can lead to close encounters that quickly destabilize the 
system. Conversely, for initial conditions close to equi¬ 
librium, conjunctions always happen at the most stable 
longitude, minimizing the strength of the kicks incurred 
at conjunctions. 

The equilibrium location is set by the planet-pair’s 
proximity to exact resonance. When the planets’ periods 
are far from resonance, the equilibrium configuration is 
for both orbits to be circular. As one brings the semima¬ 
jor axes toward the resonant separation, the equilibrium 
eccentricity shifts away from the origin toward the right 
along the x axis. This so-called forced eccentricity is 
given by 



( 11 ) 


where we have omitted a coefficient of order unity (ss 
0.8 in our case, see Eq. 11 in Goldreich and Schlichting 


90 ° 



Fig. 4. — Polar plot of e and </> for a test particle near resonance 
with an outer perturber on an exterior circular orbit. The fixed 
point’s location at ey ~ 0.008 (blue dot), is set by Eq. 11, where 
we have used fi — 10 —4 and A = 0.01 (Eq. 12). Values of (e, (f>) 
for systems with four different initial conditions trace out the red 
solid curves. As the point describing the system moves along one 
of these curves, e and <f> vary (the point lies at different azimuths 
and intersects different dashed circles of constant eccentricity). 

As long as the initial eccentricity eo is not ef (in which case the 
system traces approximately traces out a circle at constant eo), the 
eccentricity variations are comparable to ef. 

2014), /i is the ratio of the perturber’s mass to that of 
the central body, and A is a dimensionless measure of 
the proximity to resonance (Lithwick and Wu 2012) 12 





( 12 ) 


where the first quantity in parentheses approximately 
approaches zero at resonance. When planets are more 
widely separated than the resonant separation ratio, 
A > 0, and when they are closer together, A < 0. 1; 
We note that when A < 0, Eq. 11 implies e/ < 0. This is 
the correct result for one of the fixed points if one consid¬ 
ers a negative eccentricity to correspond to <j> = tt rather 
than <f> = 0. 

For rough estimates close to resonance (the region of 
interest), one can write the approximate forms 



3 + 1 

3 

(13) 

~ ° 2 

(A) . 

(14) 

ai 

\ ai J 

\ / res 


12 We have borrowed the intuitive variable A from Lithwick and 
Wu (2012), but redefined it in a way that renders Eq. 11 partic¬ 
ularly simple. Our A thus varies trivially from theirs, and agrees 
for A < 1. 

13 Equation 11 unphysically diverges at A = 0. This is due to 
a truncation of the equations of motion at leading order in the 
small quantities e and fi. This is not a limitation since, as dis¬ 
cussed below, dissipation prevents planets from approaching exact 
resonance. 























7 


where ( 02/01 ) res is the semimajor axis at nominal reso¬ 
nance [(j + l)/j] 2 / 3 . Close to resonance, the variable A is 
approximately therefore simply the shift of the period or 
semimajor axis ratio from the nominal resonance value 
(i.e., Eq. 10 with vj = 0). We further note that Eq. 14 
generalizes to the case discussed in Sec. 4.3 where the res¬ 
onant locations shift. In that case, (a 2 /ai) res would cor¬ 
respond to the semimajor axis ratio that satisfies Eq. 9. 


4.2. When are resonances important? 

We are now in a position to evaluate the effects of 
resonances on our analysis in Sec. 3.4. We found in the 
previous section that, for a given A, the fixed point in 
the dynamics lies at e/ ~ pj A (Eq. 11). Assuming the 
initial eccentricities aren’t ^ e/ (see caption to Fig. 4), 
ef sets the scale for the eccentricities the system will 
reach—orbits at the fixed point will remain at e = ef, 
while in the limiting case of initially circular orbits, the 
eccentricity will approximately vary between zero and 
2e/ (Fig. 4). 

This can explain the discrepancies in the figures of 
Sec. 3.4 from the expected mass thresholds in Eqs. 6 and 
7. Those simulations started with circular orbits. If the 
semimajor axes for a pair of planets lie within some A of 
resonance, the resultant variations in e ~ p/A become 
an appreciable fraction of the eccentricity required for 
orbits to cross e cross « A a/a. 

In our nominal 5-planet case, we expected to find the 
mass threshold at ss 3OM0, or /i = M p /M* « 2 x 10~ 4 . 
Additionally, the period ratio between the outermost two 
planets P 5 /P 4 ~ 1.32, which is close to the first-order 
4:3 MMR (with period ratio 1.33). From Eq. 12, A ss 
—0.01, so we expect eccentricities to reach e ~ 0.02 by 
Eq. 11. This is an appreciable fraction (~ 10%) of the 
eccentricity required for orbit-crossing e cross « A a/a ~ 
0 .2, which hastens the onset of instablility. 

If we move the system away from resonance by shifting 
the initial semimajor axis of the outermost planet inward 
by 1 AU, A triples, correspondingly decreasing the mag¬ 
nitude of eccentricity oscillations. Performing the same 
analysis as in Sec. 3.4 with this shifted system, we ob¬ 
tain Fig. 5, which matches well the simple expectations 
from Sec. 3. 14 This suggests an approximate empirical 
threshold for e/ of ~ 10% of the eccentricity required for 
orbit-crossing. 

To investigate this claim further, we ran a large suite 
of numerical integrations for two planets in the limits 
of a massless inner planet and of two equal-mass plan¬ 
ets. The inner planet was placed at the location of gap 
3, and the second was placed at various radii near the 
locations of the major first-order MMRs (2:1, 3:2 and 
4:3). We started the planets on circular orbits, and for 
a given proximity to a particular resonance A, we var¬ 
ied the planet mass to see at what value of e/ = p/A 
(Eq.ll) the resonantly induced eccentricity oscillations 
led to instability within 10 6 yrs. As might be expected, 
this critical e/ varies with A and the particular reso¬ 
nance; however, we find that in all cases, the critical e/ 

14 The extra stability in Fig. 5 at r e = 10 4 yrs is due to the 
damping timescale being faster than the slower instability timescale 
beyond 3.5 Rh> as discussed in Sec. 3.3. 
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Fig. 5. — Time to instability for a non-resonant case (the out¬ 
ermost planet has been shifted inward by 1 AU relative to our 
nominal setup). Otherwise same as Fig. 2. 


is less than 0.01. We can thus define (through Eq. 11) 

A crit = 100 p. (15) 

This provides a useful rule-of-thumb: For planets near 
the 2:1,3:2, and 4:3, if A > A cr it, then the resonances do 
not significantly affect stability, and Eqs. 6 and 7 should 
provide a robust mass limit for stability. This is also the 
limit in which “naive” initial conditions of circular or¬ 
bits are appropriate for a numerical stability analysis in 
an eccentricity-damping disk. Conversely, if A < A cr i t , 
stability will be more sensitive to initial conditions and a 
more detailed analytical or numerical analysis is required 
(see Sec. 5). It is important to note that this is a state¬ 
ment about stability over short timescales of 10 3 — 10 4 
orbits, and that, particularly for the 4:3, one should en¬ 
sure that the A values one considers does not place the 
system closer to a different first-order MMR. 

As a second example, in our four-planet case (Fig. 3), 
the expected mass was approximately half a Jupiter 
mass, or /i ss ICR 3 , and A « 0.005. We therefore ex¬ 
pect eccentricitiy variations ~ 0.1, which should be very 
unstable, and helps explain the large deviations from 
the expected thresholds in Fig. 3. From Eq. 15 we re¬ 
quire A > A cl .jt = 0.1. This can be achieved by shift¬ 
ing the second-to-last planet outward by 4 AU. Figure 6 
shows the result, which matches very closely the analytic 
threshold from Eq. 6. One should compare the shifts re¬ 
quired to move between resonant and non-resonant cases 
to the observational uncertainties in the orbital radii. 
But before doing so, we must consider an important com¬ 
plication. 

4.3. The resonances are not where you think they are 

The locations of resonances in a system are usually 
accurately determined by simple period ratios, and are 
thus independent of both the central star’s mass and 
the absolutely size scale of the system. However, this 
is only true when pericenter precession rates are slow 
compared to the orbital rates (see Eq. 10). This is gen¬ 
erally a safe assumption, as the pericenter direction is a 
conserved quantity in the Kepler problem, so small non- 
keplerian perturbations Enk induce a slow precession of 
order w/n ~ Fnk/Pk "C 1, where Fk is the dominant 
gravitational force from the central star. 
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Fig. 6.— Time to instability for a non-resonant case of our 
4-planet system (with the planet with co-orbital material moved 
from 71.2 to 75.5 AU). 

However, for young systems, the disk’s gravity can 
present a significant non-Keplerian component. In the 
case of HL Tau, Kwon et al. (2011) estimate a disk mass 
Md of 0.13 M 0 , which represents about one quarter of 
the central star’s mass. Assuming a razor-thin, axisym- 
metric disk, this would induce a pericenter precession 
rate (relative to the mean motion) of w/n ~ Mp /M* ~ 
0.25 (Eq. Bll). According to Eq. 10, this would move 
the location of the 4:3 resonance by ~ 10%. This is a 
significant shift, since we found in the previous section 
that in our nominal five-planet model, the 4:3 resonance 
between the outermost two planets was only important 
within A ~ 1% of resonance. 

We point out, however, that the precession rate derived 
in Eq. Bll is very uncertain, as it ignores deviations from 
a smooth axisymmetric disk. These can be particularly 
large in the vicinity of the planet (e.g., any gap or spiral 
density waves induced by the body), and would provide 
contributions to the precession rate that are difficult to 
calculate due to the feedback between the planet and 
disk, the dependence on disk structure and thermody¬ 
namics, etc. Nevertheless, Eq. Bll provides the scale on 
which the precession rates (and therefore the resonance 
locations) are uncertain. We also note that this base 
pericenter precession rate induced by a smooth axisym¬ 
metric disk is negative, which means through Eq. 10 that 
resonances lie at wider period ratios, and thus at wider 
separations, than one would expect in the disk-free case. 
This means that the outermost two pairs of planets in our 
nominal five-planet case could plausibly lie in 5:4 reso¬ 
nances, despite their naive period ratios (using Kepler’s 
3rd law) putting them closer to the 4:3. 

We note that in contrast to the simple case consid¬ 
ered in Sec. 4.1 where the inner planet is massless, in the 
general case of two massive bodies, both orbits evolve, 
and there exists a second resonant argument involving 
the longitude of pericenter of the outer body. Beauge 
et al. (2003) (see also Beauge et al. 2006) find that in 
the disk-less case, the equilibrium solutions in a secular 
approximation correspond to the two percenters process¬ 
ing at the same rate (apsidal corotation). An interesting 
question we leave for future work is whether in a massive 
disk, the differential pericenter precession rates induced 
by the disk are strong enough to break the apsidal corota¬ 
tion that is expected from slow resonant capture (Beauge 


et al. 2006) . Additionally, this interaction between plan¬ 
ets acting in addition to the disk should affect the precise 
location of the shifted resonances; however, because the 
disk-induced precession is so uncertain, we do not pursue 
this further. 


4.4. Can one observationally constrain whether ALMA 
planets are in resonance? 

The previous section argued that the locations of res¬ 
onances are uncertain for planets in a disk by of order 
the ratio of the disk’s mass to the central star. While 
this ratio is typically observed to be ~ 1% (see Fig. 5 in 
Williams and Cieza 2011), there is a large scatter of ~ ±1 
dex, and one might expect early ALMA observations to 
target brighter disks. This suggests that in many sys¬ 
tems, even if one could extract exact semimajor axes for 
each of the planets, one could not determine the proxim¬ 
ity to resonance A to within the accuracy needed to tell 
whether the resonant interactions are strong (Eq. 15). 

Additionally, there are observational uncertainties in 
the orbital radii due to the spatial resolution. From 
Eq. 14, the error in A (6 A) associated with an error Sa 
on the semimajor axis is 



(16) 


In the HL Tau image, the resolution is about 5 AU 15 , 
which corresponds to an uncertainty in A of sa 0.1. 

Put together, this suggests that one is unlikely to deter¬ 
mine whether planets in ALMA systems are resonantly 
interacting from observed gap radii. Conversely, it will be 
similarly difficult to rule out resonances in tightly-packed 
systems (all semimajor axis ratios less than 1.4 are within 
ss 10% of a first-order MMR). However, a promising av¬ 
enue for observationally breaking this degeneracy is to 
search for eccentric gaps (see Sec. 5.4). 


5. UPPER MASS LIMITS IN ALMA SYSTEMS 

Given the above uncertainties, we proceed to estimate 
the maximum masses across the range of possibilities. 
This will correspond to systems that are stabilized near 
the centers of resonances (Fig. 4). We have already seen 
that in resonant configurations, different initial condi¬ 
tions can lead to very different outcomes. This presents 
a large parameter space to consider; however, dissipa¬ 
tion from interactions with the disk powerfully shrinks 
the relevant phase space. 


5.1. Interactions with the disk 

Goldreich and Schlichting (2014) derive an analytic 
theory for the parametrized disk-damping model used 
in this paper, assuming the damping timescales are long 
compared to the resonant timescales. This is only strictly 
satisfied (at our lowest masses) in the simulations with 
r e = 10 6 yrs, but their analytical formulae nevertheless 
provide a useful lens through which to interpret our re¬ 
sults. 

Eccentricity damping acts to damp the system, not 
necessarily to e = 0, but rather to the nearest fixed point 

15 https://public.nrao.edu/news/pressreleases/ 

planet-formation-alma 
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in phase space. If planets are far from resonance, so 
e/ « 0 (Eq. 11), then this indeed damps systems onto 
circular orbits. On the other hand, near a resonance, 
eccentricity damping will move the system to the forced 
eccentricity, i.e., the most stable configuration. 

At the same time, disks that damp eccentricities should 
also induce radial migration. We implement a simple mi¬ 
gration scheme by adding a force in our N-body simula¬ 
tions counter to the particle’s velocity that continually 
saps orbital energy (Papaloizou and Larwood 2000), 



where F a is the semimajor-axis-damping force vector act¬ 
ing on a planet, and v is the body’s velocity vector. 
This yields an orbit-averaged exponential decay of the 
semimajor axis on a timescale r a . Following Lee and 
Peale (2002), we impose a simple proportionality rela¬ 
tion r a = Kr e , with a nominal value of K = 100. This 
is consistent with the theory for planets embedded in a 
disk (Artymowicz 1993; Ward 1997) (Type I migration), 
and may apply to planets that have opened a partial 
gap. For simplicity, our simulations apply the migration 
force only to the outer planet, so r a shoud be viewed as 
the timescale for relative migration, which is the relevant 
quantity for the resonant evolution (e.g., Goldreich and 
Schlichting 2014). In order for differential migration to 
lead to capture in resonance, planets must start wide of 
resonance (A > 0) and migrate convergently, rather than 
vice-versa (e.g., Murray and Dermott 1999, Sec. 8.12 and 
8.13). In such a case, convergent migration pushes plan¬ 
ets deeper into resonance on a timescale r a , i.e., it moves 
the fixed point to the right in Fig. 4 (see Eq. 11). 

This would likely eventually lead to instability, but ec¬ 
centricity damping provides a balancing effect, by push¬ 
ing the stable conjunction longitude away from </> = 0. 
This is in some ways analogous to how dissipation in the 
Earth, coupled with its rotation, causes the tidal bulge 
raised by the Moon to lead the axes connecting the two 
bodies. In both cases, this asymmetry allows the objects 
to secularly exchange angular momentum, and causes the 
two bodies to separate from one another. This “reso¬ 
nant repulsion” has been quantitatively investigated by 
Greenberg (1981) and more recently by Lithwick and Wu 
(2012) and Batygin and Morbidclli (2013) to explain the 
overabundance of planet pairs wide of resonance in the 
Kepler mission’s sample. This orbital divergence causes 
A to increase, which would move the fixed point toward 
the origin in Fig. 4, following Eq. 11. Because resonant 
repulsion is stronger closer to resonance, the system ap¬ 
proaches an equilibrium forced eccentricity e eq where the 
migration is balanced by resonant repulsion. This occurs 
at 

e eq ~ = A'" 1 / 2 , (18) 

where in general, this equlibrium eccentricity will depend 
not only on the assumed K, but also on the mass ratio 
between planets, and the balance of energy and angu¬ 
lar momentum dissipation employed for the eccentricity 
damping (see Papaloizou and Szuszkiewicz 2005; Goldre¬ 
ich and Schlichting 2014). Regardless of the exact val¬ 
ues, dissipation shrinks an otherwise large phase space of 


initial conditions to tractable equilibrium configurations 
amenable to numerical investigation. 

In such systems acted on by migration and eccentricity 
damping, Goldreich and Schlichting (2014) further find 
that as one decreases the perturbing planet’s mass be¬ 
yond a critical threshold, the fixed point changes stability 
and pushes the system away from the forced eccentric¬ 
ity. In their approximate analytical model, this can cause 
the system to lose the resonance, in which case the plan¬ 
ets migrate through the resonance configuration and do 
not capture into the equilibrium configuration discussed 
above. We do not see such outcomes in our simulations. 
This may be due to our working outside their regime of 
approximation, i.e., we integrate equal-mass planets, and 
in some cases the damping timescales are comparable to 
the resonant (libration) timescales. 

5.2. Stochastic Migration 

The discussion above assumes a smooth proto¬ 
planetary disk which leads to smooth migration and ec¬ 
centricity damping forces. Real disks are expected to be 
turbulent (Nelson 2005) either due to the magneto rota¬ 
tional instability (MRI) or the gravitational instability, 
and this can cause planets to migrate (e.g., Nelson 2005; 
Johnson et al. 2006; Rein and Papaloizou 2009). 

The precise level of turbulence is still subject to debate 
and might depend on a variety of factors such as surface 
density, surface density gradient, cooling rate, local ion¬ 
ization fraction and dust properties, but an exhaustive 
discussion of these effects go beyond the scope of this 
paper. Here, we parametrize the effect of a non-uniform 
surface density due to any kind of turbulence with the 
magnitude of density fluctuations in the disk only: 

= AE/E. (19) 

In numerical simulations, this value has been evaluated 
to be « 0.1 in the case of MRI turbulence, in the case 
of dead zones it is weaker, perhaps of order 0.01, with 
significant uncertainties remaining (Nelson 2005; Oishi 
et al. 2007; Nelson and Gressel 2010; Gressel et al. 2011; 
Yang et al. 2012). A value closer to unity can be expected 
from strongly gravitationally unstable disks (Rice et al. 
2011). 

The specific force exerted by a density perturbation 16 
on a nearby planet can then be written as 

Aas ~ -ttGS^T, (20) 

For many types of instabilities in the disk, absent long- 
lived coherent structures such as vortices, one might ex¬ 
pect the correlation time of density perturbations to be 
comparable to the orbital timescale, r c ss n -1 . Hav¬ 
ing the magnitude of random forces and their correlation 
time, we can define a diffusion coefficient D to completely 
describe the effects of random forces: 

D = 2 F 2 s t c = U 2 G 2 <5 2 Y?n-\ (21) 

Using the results of Rein and Papaloizou (2009), we can 
now estimate the distance A a that a planet migrates in 

16 Note that this expression is independent of the perturbation 
size ax; since the mass of the perturbation scales as ~ 
and the force on a planet scales as GM ( 5 x;/a|,. 
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a time t purely due to stochastic forces: 

A a = yfADt/n 2 (22) 

Note that the \Jt dependence implies that we are model¬ 
ing this process as a simple random walk in a. Assuming 
a power-law dependence with distance, £ oc r -7 (7 < 2 ), 
with a sharp cutoff at r = R Cl and relating its normal¬ 
ization to the mass of the disk (see Eq. B3) , we can write 
down the relative change in a in a compact form as a 
function of time, disk mass and magnitude of density 
fluctuations: 



Planets close in are affected more than planets further 
out due to their shorter dynamical timescale. Assuming 
the value for 7 traditionally assumed for the minimum- 
mass solar nebula of 1.5 (Hayashi 1981), and using a 
conservative value of 5s =0.01, an HL Tau disk mass of 
Md/M* ss 0.25 (Kwon et al. 2011) and a period for the 
innermost planet of ~ 100 yrs, we find that within one 
million yrs, the planet can move by A a/a « 0.1. 

This result suggests that stochastic forces from a tur¬ 
bulent disk as massive as that of HL Tau could move 
around planets significantly during their lifetime. How¬ 
ever, stochastic migration could be mitigated if planets 
open deep gaps in the gas disk (Rein and Papaloizou 
2009) and, depending on the putative planets’ masses, 
other migration mechanisms (e.g., Type I Ward 1997) 
could dominate the migration. In any case, with ALMA 
and other observatories putting limits on several of the 
factors in Eq. 23, it will be an exciting task to constrain 
the effects of turbulent migration directly from disk ob¬ 
servations (see, e.g., Hughes et al. 2011; Simon et al. 
2011; Shi and Chiang 2014; Simon et al. 2015). 

5.3. Initial conditions 

If planets migrate smoothly toward a particular reso¬ 
nance, the system should remain near the fixed point as 
the forced eccentricity rises from « 0 far from the res¬ 
onance to the equilibrium eccentricity. In this scenario, 
when planets are near resonance, initially circular or¬ 
bits are therefore not the appropriate initial conditions. 
They unphysically result in large-amplitude oscillations 
in both eccentricity and conjunction direction </> (Fig. 4). 
This explains why resonances lowered our mass limits in 
Figs. 2 and 3, which used initially circular orbits. 

Conversely, initial conditions at the fixed point should 
promote stability and raise our mass limits. While we 
have so far considered a simplified dynamical model to 
gain intuition, in the real problem, each orbit will have 
its own forced eccentricity that will depend on the differ¬ 
ent planet masses, damping timescales and proximities 
to various resonances. It is therefore impractical to ap¬ 
propriately assign initial conditions. Instead, we allow 
the eccentricity damping to move the system to the fixed 
point by numerically migrating the planets into equilib¬ 
rium configurations ahead of the “start” of our simula¬ 
tions. 

Of course, if the planetary masses are large enough that 
the separations correspond to less than the 3.5 Rh limit 
(see Fig. 2 and Eq. 6 ), the instability time is short, so 


before planets can migrate into the stable resonant con¬ 
figuration, the system will disperse. However, the plan¬ 
ets are growing throughout this process. It is therefore 
plausible for planets to migrate into resonance at lower 
masses (when their separations correspond to more than 
3.5 Rh) and grow together in this stable configuration 
beyond the limit imposed by Eq. 6 . This should provide 
a robust upper limit to the masses. 

This effect is illustrated in Fig. 7, which shows a sim¬ 
ulation of cores growing from 10 M® (time increases up¬ 
ward in the plot). The system of planets in blue has been 
migrated into a chain of 4:3 resonances among the outer 
three bodies, while the system of red planets has not. 
Once the planets grow above the dotted horizontal line 
of 45 M® (the 3.5 Rh boundary for these separations), 
the non-resonant planets (red) destabilize on the short 
timescales on which conjunctions occur (Sec. 6 ), while 
the blue planets continue to grow, protected from close 
encounters by the resonance. Eventually the masses grow 
to the point that kicks at conjunctions knock the planets 
out of resonance, and swift instability ensues. 



semimajor axis [AU] 

Fig. 7.— Two systems of planets, growing from 10 cores 
(time increases upward). The outer three blue planets are in a chain 
of 4:3 resonances, while the red planets are not. The horizontal 
dotted line represents the mass threshold above which the orbital 
separations represent less than 3.5 Rh-> and non-resonant systems 
(red) should rapidly destabilitize. The resonant (blue) system can 
grow to larger masses through the resonance mitigating the effects 
of close encounters. 

We therefore begin with giant planet cores of 10M®, 
and migrate them into the closest resonant configura¬ 
tion. 1 ' In our nominal 5-planet case, both of the out¬ 
ermost pair of planets are near resonance, with period 
ratios of 1.29 and 1.32, respectively. We therefore start 
them at wider separations 18 , and apply inward migra¬ 
tion on the outermost planet at a rate r a = 100 r e until 
they capture in a chain of 4:3 resonances. Once the re¬ 
spective resonant angles librate with small amplitude, we 
adiabatically (on a timescale 10 times longer than the li- 
bration period) raise the masses to the value in question. 
From this point, we then measure the time to instability 
as before. 

1 ' In the r e = 10 4 case, the migration time is fast enough that we 
have to start with M = 40M® in order to capture into resonance. 
See Ogihara and Kobayashi (2013) for how the mass required for 
capture is related to the migration timescale. 

18 We numerically selected our semimajor axes to yield the ob¬ 
served orbital radii at the end of this “warmup” of initial conditions 
(to within 2.5 AU, or half the image’s spatial resolution). 
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One should ask whether this warm-up time is signif¬ 
icant. We are effectively assuming that planets start 
in the resonance fully formed, but planets take time to 
grow, and the equilibrium configuration takes time to set 
up. How much time this represents depends on uncertain 
paramaters such as how close the planets start initially, 
how the migration timescales vary as the planet grows in 
mass, etc. A coherent picture for a system like HL Tau 
would have to consider these issues. For simplicity, we 
sidestep this complication, effectively assuming that the 
time for planets to grow and for orbits to migrate into res¬ 
onance is much shorter than the age of the system. This 
means that our estimates provide a conservative lower 
limit for how large planets can grow while in resonance. 

5 . 4 . Maximum masses in resonant configurations 

The results for our nominal five-planet case are shown 
in Fig. 8. Areas in white captured into resonance and 
were stable for 1 Myr, while areas in black were too un¬ 
stable to even reach an equilibrium configuration. The 
transition is abrupt, and there is only a slight depen¬ 
dence on r e , which suggests that we can robustly set a 
mass threshold of M ss 95M®, or about one Saturn- 
mass. We reiterate that this represents the maximum 
mass that can survive for at least 1 Myr in a resonant 
configuration. In reality, planets take time to grow and 
migrate into resonance, so this implies the planet masses 
could be somewhat greater than our estimate (how much 
greater depends on uncertain parameters). This can be 
seen in Fig. 7, where the planets grow « 40% beyond a 
Saturn mass before destabilizing. We also point out that 
the resonant configuration pushed the mass limit signif¬ 
icantly beyond the 3.5 Rh criterion (dashed blue line in 
Fig. 8), see Eq.6). 
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Fig. 8.— Maximum masses for planets in the HL Tau system. 
This assumes a scenario where planets grow while in equilibrium 
resonance configurations (see Sec. 5.3 for details). 

Performing the same exercise on our more widely sep¬ 
arated 4-planet case, we find that planets can grow to¬ 
gether in resonance to at least 230 M®. 

An important consequence of resonant interactions is 
that they force eccentricities that will not be damped 
away. If gaps in nascent systems are observed to be ec¬ 
centric, and the period ratios are consistent with a res¬ 
onant configuration (Sec. 4.4), that would provide tan¬ 
talizing evidence for the interpretation that planets are 


carving the gaps, and would break the degeneracy as to 
whether planets are resonant or not. 

The magnitude of the equilibrium eccentricities in 
our assumed model will depend on both K (Sec. 5.1) 
and the mass ratios between planets (Papaloizou and 
Szuszkiewicz 2005; Goldreich and Schlichting 2014). In 
our simulations with K = 100 and equal masses, we find 
the equilibrium eccentricities for the gap 3 and gap 5 
planets are approximately 0.015, and e « 0.03 for the 
gap 4 planet, independent of r e and the assumed mass. 
A survey of the parameter space is beyond the scope 
of this study, but observed eccentric gaps would clearly 
motivate such a study. Smaller values of the uncertain 
parameter K would yield larger equilibrium eccentricities 
(see Eq. 18). 

The above discussion implicitly assumes that a planet’s 
orbital eccentricity is communicated to the gap it carves. 
One might expect that, as long as the gap carving time 
is short compared to the timescale on which the planet’s 
orbit changes, this should approximately hold true. This 
is certainly an important question meriting future study. 

An additional complicating factor in the case of HL 
Tau is that the massive disk should cause the pericenters 
to precess on timescales (~ a few orbits, see Sec. 4.3) 
that are likely comparable or faster than the timescale on 
which planets carve gaps (~ 10 — 100 orbits Fouchet et al. 
2010). This would smear out any eccentric signature into 
a circular gap. This provides a strong motivation for 
probing more typical, low-mass disks, Mp/ A/* < 1%, 
where the precession periods are likely longer than the 
timescale for gap opening, and eccentric features could 
be preserved. 


6. CONCLUSIONS 

We divide our conclusions in two. First we summarize 
how one might analyze the dynamical stability of an ob¬ 
served system, and then we recapitulate our findings for 
HL Tau, assuming the gaps are carved by planets. 

6 . 1 . Dynamically constraining ALMA systems 

In this paper, we have provided a framework for un¬ 
derstanding the stability of general planetary systems in¬ 
teracting with a gas disk. The short disk lifetimes (< 3 
Myr) over which one might detect nascent systems al¬ 
low for particularly simple estimates. 19 We now outline 
a simplified workflow for analyzing a new system, and 
reference the more detailed discussion in the body of the 
paper. These dynamical constraints can easily be incor¬ 
porated to guide the design and interpretations of both 
observers and those simulating gap opening in observed 
systems. 

Our treatment has assumed equal-mass planets. In re¬ 
ality, our procedure effectively estimates the maximum 
masses M cr ; t for the most closely spaced planets in the 
system (in terms of A a/a). In principle, a similar analy¬ 
sis carried out for the progressively more separated pairs 
of planets in the system should yield reliable results. It 
would be an interesting extension of this work to test how 

19 The maximum masses so obtained also apply over longer 
timescales, they just become less informative. On Gyr timescales, 
secular interactions (e.g., Lithwick and Wu 2011; Wu and Lithwick 
2011; Batygin et al. 2015) and slower chaotic diffusion can set much 
more stringent mass limits, at the cost of considerable complexity. 
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accurately these simple criteria can be applied in this it¬ 
erative fashion. We now summarize a simple procedure 
for analyzing a system’s stability. 

• The base estimate for M cr j t is given by inserting the 
minimum fractional separation among planet pairs 
(Aa/a) m i n in Eq. 6. If these planets are not near 
resonances (see below), this provides a stringent 
upper mass limit for the planets in question, which 
is independent of the assumed interactions with the 
disk. 

One can then proceed to refine this limit, which only 
rigorously applies in the case of a two-planet system. 

• If planets on either side of the closest two have 
A a/a < 4(Aa/a) m in, and the eccentricity damping 
timescale is long compared to the time betwen con¬ 
junctions (see footnote in Sec. 6 ), planet masses be¬ 
low the above-calculated will destabilize over 
the system’s lifetime (Sec. 3.2). 

• In this large r e regime, and in the limiting case that 
adjacent pairs share A a/a = (Aa/a) m i n , M cr ; t is 
pushed down to the value given by Eq. 7. In inter¬ 
mediate cases, Eqs. 6 and 7 bracket M cr j t (Cham¬ 
bers et al. 1996; Chatterjee et al. 2008; Marzari 
2014). 

Near-resonant orbits will also modify M cr j t . 

• The locations of resonances are uncertain to within 
~ Mdi s k/M* (Sec. 4.3), and the orbital radii to 
within some observational uncertainty (Sec. 4.4). 

• If, within these uncertainties, the orbits’ fractional 
distances from first-order MMRs are always large 
compared to M crit /M* (Sec. 4.2), then M crit is a 
robust estimate. 

• If first-order MMRs can not be ruled out, one can 
numerically set an upper mass limit by assuming 
the planets have reached an equilibrium resonant 
configuration (Sec. 5.3) 

6 .2. HL Tau 

Our main conclusions are summarized in the bullet 
points below , 20 

• A disk’s mass should shift the resonance locations 
from integer period ratios by ~ M^/M* (~ 10 % 
for HL Tau), but the precise value is uncertain due 
to theoretical difficulties in evaluating the contribu¬ 
tion to orbital precession rates from local material 
disturbed by the planet (Sec. 4.3). 

20 Note added during revision: An analysis of the HL Tau data 
has recently appeared (ALMA Partnership et al. 2015). They find 
that gaps 3 and 4 are closer than we assume; using their values 
would somewhat lower our quoted maximum masses in the 5-planet 
case. ALMA Partnership et al. (2015) also use the disk kinematics 
to derive a mass that is approximately twice as large (~ 1.3M@) as 
inferred by previous studies. If indeed the star is this massive, and 
the disk is somewhat less massive than suggested by the modeling 
of Kwon et al. (2011), this would correspondingly reduce the shifts 
of resonances from their nominal locations (Sec. 4.3), and mitigate 
the rapid erasure of eccentric gaps by precession (Sec. 5.4). Re¬ 
markably, ALMA Partnership et al. (2015) find each of the gaps 
to be offset from center. It will be an exciting task to determine 
whether these offsets are due to eccentric gaps and to exploit this 
information to further probe the planetary hypothesis. 


• This theoretical uncertainty, coupled with observa¬ 
tional uncertainties in the orbital radii make it dif¬ 
ficult to constrain whether planets are resonantly 
interacting. Conversely, in closely spaced systems 
like HL Tau, one cannot rule out resonant configu¬ 
rations (Sec. 4.4). 

• Despite the possible window into the era of planet 
formation afforded by the HL Tau data, these un¬ 
certainties in the dynamical state are an important 
consideration when using the HL Tau system (or 
others discovered in the future) as initial conditions 
for simulations of, e.g., planet scattering. 

• Observationally detecting eccentric gaps that may 
be resonant would provide tantalizing evidence 
that resonantly interacting planets are causing 
them. This would additionally break the degen¬ 
eracy in the system’s dynamical state, constrain 
the pericenter precession rates and provide proxies 
for initial conditions at the end of planet formation 
(Sec. 5.4). 

• The rapid pericenter precession rates induced by 
massive disks will erase eccentric signatures and 
generate circular gaps (Sec. 5.4). This strongly mo¬ 
tivates pushing toward less massive, more typical 
disks with Mp/M* < 10 -2 . 

• The timescale for stochastic diffusion onto crossing 
orbits, as well as observational limits on the level 
of density fluctuations in the disk, can set limits on 
the level of turbulence in the outer disk regions of 
systems probed by ALMA (Sec. 5.2). 

• Planets can grow to much larger masses than one 
would naively expect if they capture in resonance 
at low masses and grow together (Fig. 7). This can 
continue until the masses become large enough that 
kicks at conjunctions knock the planets out of reso¬ 
nance, leading to swift instability. It would be valu¬ 
able to investigate how the threshold mass varies 
between different first-order resonances, as well as 
with the assumed dissipation timescales. If dis¬ 
tant planets tend to capture into resonances, this 
could represent a universal pathway to instability 
(Sec. 5.3). 

As an illustration of the framework we have elaborated, 
we focused primarily on a nominal model with five plan¬ 
ets orbiting in the most prominent gaps (Fig. 1), as well 
as a 4-planet model where gaps 3 and 4 are generated by 
a single planet lying between the gaps, with co-orbital 
material along its path. In both cases, the inner two 
planets are largely dynamically decoupled from the sys¬ 
tem, so the below mass limits apply to the outer planets. 

• 5-planet Case (Fig. 1): If planets are not resonantly 
interacting, the maximum masses that allow for 
stability over the system’s age are « 2 Neptune 
masses. If planets grow together while in reso¬ 
nance, masses can reach at least 1 Saturn mass. 

• 4-planet Case: In the non-resonant case, planet 
masses can be as large as ~ 1 Saturn mass, while 
growing together in resonance, masses can reach at 
least 230 M e . 
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• The inner planets are largely decoupled from the 
system and are consistent (dynamically) with at 
least several Jupiter masses (Eq. 7). At these large 
masses, higher-order corrections in the planet-to- 
star mass ratio become important, and our equa¬ 
tions may not be precise. 

Additional observations will likely produce high resolu¬ 
tion maps of the gas disk. Whether gas is observed at the 
locations of gaps in the dust will additionally constrain 
the maximum masses of any planets carving them. In 
addition, comparing the dust and gas distributions can 
help pinpoint the stage of planetary formation the system 
is in: during the build-up of giant-planet cores, during 
the runaway accretion of gas, or in between both phases. 


Given constraints on the system’s age, it may soon be 
possible to empirically test theoretical expectations (e.g., 
Pollack et al. 1996; Alibert et al. 2005). Putting new ob¬ 
servations and simulations in the dynamical context we 
have presented should prove an exciting and fruitful task. 

We would like to thank Yanqin Wu, Norman Murray, 
Hilke Schlichting, Pawel Artymowicz, Ruth Murray-Clay 
and Katherine Krettke for insightful discussions. We are 
also grateful to the anonymous referee who helped im¬ 
prove this manuscript. This research was made possible 
by the Sunnyvale cluster at the Canadian Institute for 
Theoretical Astrophysics. 


APPENDIX 

APPENDIX A: DETERMINATION OF ORBITAL RADII FROM THE ALMA IMAGE 

Kwon et al. (2011) used CARMA 21 observations in a band centered at the same wavelength (1.3 mm) as the 
ALMA image to estimate the disk’s diameter at 235 AU. Additionally, they constrained the PA angle at 136° and the 
inclination of the disc at 40°. Using commercial software, we drew ellipses onto the press-release image. We adopted 
a PA angle of 136°; however, through trial and error, we found that if one assumes the orbits are circular, the aspect 
ratio between the minor and major axes implies an inclination of the system of « 48°. This is likely caused by a non 
circular beam on the ALMA image, which we could not compensate for. The minor axis only served the purpose of 
centering the ellipses and does not impact the value that interest us: the orbital diameter of each of the gaps, which 
we estimated using the major axis. 

We started by forcing all ellipses to share a common center; the aspect ratio was also kept constant. The width of 
each ellipse’s circumference was set to correspond to the average width of the gaps. This meant that the inner and 
outer edges both contribute to identifying the middle of the gap. We visually adjusted only the major axes, each 
independently. A first ellipse was adjusted to one of the gaps (randomly chosen). The order with which subsequent 
gaps were located, was random. The full procedure was repeated five times. This means that each gap has five orbital 
diameter estimates. We gathered the set and computed the mean of each gap’s diameter and provide an error using 
the rms of the values. 

In Fig. 1, gap 5 is visibly not well centered. Because of this, a second series of five adjustments was made, reproducing 
the previous methodology except that now the center of each ellipse was independently positioned. We obtained 
consistent results within the error bars, so we conclude that any eccentricity in the outer ring is not significant at the 
precision level of our approach. 

Both sets being in agreement, each gap has ten estimates, which we combined to produce a single set of orbital 
distances. We thus obtain diameters in arbitrary units, which we denote pt. To convert to an absolute scale using 
the CARMA diameter of 235 AU (Kwon et al. 2011), we estimated the diameter of the disk by drawing an additional, 
thin, ellipse, ten times independently. This yielded a diameter of the disk of 1385 ± llpt, which we used to convert 
the orbital radii into AU in the second column of table A. 

We do not include in our estimates an error incurred by this unit conversion (the absolute scale is uncertain to at 
least the level of the error in the distance to HL Tau, approximately 15% Torres et al. 2007), partly because this 
uncertainty in absolute scale should not have a strong impact on our results. The major effect is to change the orbital 
periods for a given assumed mass for the central star. But because instability timescales are such a steep function 
of the relative separations (Sec. 7), a slight change in the number of orbits executed does not qualitatively influence 
the results. The important quantities in our analysis are the relative separations, which are well constrained by our 
method. 


TABLE 1 

Gap locations 3, as estimated from Figure 1. 

3 Errors do not include uncertainties in the image’s overall scale. Such a rescaling would have little effect on our results (see Appendix A). 


gap 

number 


diameters 
(arbitrary unit) 


physical 

orbital radii (AU) 


160.1 ± 2.6 
392.8 ± 2.2 
767.5 ± 7.5 

911.3 ±4.4 
1096 ± 11 


13.6 ±0.2 

33.3 ±0.2 
65.1 ± 0.6 

77.3 ± 0.4 
93.0 ±0.9 


21 Combined Array for Research in Millimeter-wave Astronomy 
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APPENDIX B: PRECESSION INDUCED BY A THIN AXISYMMETRIC DISK 


We consider a razor-thin, axisymmetric disk with a power-law surface density £(r), 


£(r) = £ 0 


R r 


(Bl) 


where r is the radial distance, £o is a normalization, R c is a characteristic radius, and 7 is the power-law index. If the 
disk is infinite, one obtains the simple result for the radial force Fjj(r) (Lemos et al. 1991; Evans and Read 1998), 


F D (r) = -F 0 l ^ 


(B2) 


where F 0 is another normalization. Such a disk has infinite mass, so we generalize the above solution to a case where 
one imposes a sharp cutoff at R c . Specializing to 7 < 2 (one could also add an inner cutoff to avoid divergence at 
r = 0 for 7 > 2), the mass of the disk Mjj is 


Md = 


R, 


2nY, 0 R 2 


— I (2 irr)dr = 2 _^ 


(B3) 


In this case, as in the infinite case, one can elegantly solve the problem by constructing the surface density from thin 
homoeoids 22 that are completely flattenened into a 2-D disk. The result is (see Eq. 2.157 in Binney and Tremaine 
2008), 


F D {r) 


m 


4 G f r , 

= -Jo d “ 




yV 2 — a 2 da J a y/r' 2 — 


(B4) 


where G is the gravitational constant and a represents a homoeoid’s semimajor axis. If the surface density cuts off at 
r = R Cl R c replaces the upper limit in the second integral, yielding the correction to Eq. B2 


F D (r) = -F 0 ( ^ j [l + r?(r,7)], 


(B5) 


where 


F 0 =€ 


GM d 
R 2 ' 


S = (2 - 7 ) 


n^]m ' 


T is the Gamma function, 


r l( r > 7 ) = 


2-7 (Rc 

7r( 1 r 


1-7 


2 K\ 


m 


-r br( 7 Y) 


and FI is the regularized, generalized hypergeometric function 

3 hJ{ 1/2,1/2, 2=1}, {1 ,^}, 



T[l±2] 


(B6) 

(B7) 

(B8) 


(B9) 


From Eq. B5, ry(r), represents the fractional correction due to the disk truncation, relative to the case of an infinite 
disk (Eq. B2), and rj(r) is always greater than zero. This reflects the fact that the outer material that would have 
contributed a positive force is now missing. This correction increases as r approaches R c , but is normally small. For 
7 = 1.5, the traditional value for the minimum-mass solar nebula (MMSN, Hayashi 1981), and for an HL Tau value of 
R c = 117.5 AU, 77 ~ 10% at the outermost gap. 

We can now obtain the pericenter precession rate through 


n7 

n 


Fo 1 dFu /dr 

-1 —r - 

F k 2 F k 


(BIO) 


22 The generalization of spherical shells for oblate spheroids. It 
can be shown (see, e.g., Binney and Tremaine 2008) that a mass 


lying interior to a homoeoid feels no force (a generalization of New¬ 
ton’s famous result). 
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where Fk 


GM+/r 2 is the Keplerian force from the star. This yields 


ZD 


n 


M D 

M* 


(1 



7-2 


1 + ?72(y,7) , 


(Bll) 


where 772 7) is the fractional correction relative to the case of an infinite disk (in which case one should replace Mp 

according to Eq. B3 and R c becomes a characteristic, rather than a cutoff radius). Eq. Bll gives the same functional 
dependence found by Rafikov (2013), who specialized to the classical Mestel disk, i.e., 7 = 1. 

Rather than provide an inscrutable expression for 772, we give a qualitative analysis. We first note that the two 
quantities in the square brackets of Eq. BIO have opposite signs; since both Fp and Fk are negative, the first term is 
positive, while, because Fp approaches zero at large r, dFp/dr > 0 and the second term is negative. We also point 
out that the final answer for the precession rate (Eq. Bll) is negative (recall that we have specialized to 7 < 2). If 
we now consider the effect of truncating the disk, as mentioned above, this makes Fp more negative, and the effect 
increases as one approaches R c (i.e., dFp/dr decreases, since Fp rises slower). Thus, the first term in square brackets 
in Eq. BIO becomes more positive, and the second less negative, so truncating the disk always acts to enhance the 
magnitude of zd. 

In this case, the fractional correction 772 can be quite large, due to the strong effect on the derivative in Eq. BIO. 
Again for 7 = 1.5 and R c = 117.5 AU, 772 ~ 90% at the outermost gap. Thus, while this is an important quantitative 
correction, one nevertheless obtains qualitatively correct lower limit to the precession rate by taking 772 = 0, i.e., by 
assuming an infinite disk for the force law (Eq. B2), but still using the mass of the disk for Mp in Eq. Bll. 
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